#Replication data for "The Promises and Pitfalls of 311 Data", Urban Affairs Review
#Ariel White and Kris-Stella Trump
#email Ariel (arwhi@mit.edu) with questions
#October 2016

rm(list=ls())
library(foreign)
library(data.table)

#setwd("/nfs/projects/c/civic_sp/nyc311/311asoutcome/Replication codes 311/Data preparation and merge code") #set your own WD here

NYCdata1 <- read.csv("merged_NYCcensustractdata_poldonations.csv")
NYCdata2 <- read.csv("merged_NYCcensustractdata_poldonations_publiconly.csv")
NYCdata3 <- read.csv("merged_NYCcensustractdata_poldonations_streetonly.csv")
names(NYCdata1)

NYC1 <- NYCdata1
NYC2 <- NYCdata2
NYC3 <- NYCdata3

#Some variable recoding
names(NYC1)
#Population over 18 in 2010 Census
NYC1$over18 <- NYC1$Tot_Population_CEN_2010 - NYC1$Pop_under_5_CEN_2010 - NYC1$Pop_5_17_CEN_2010;summary(NYC1$over18)
NYC2$over18 <- NYC2$Tot_Population_CEN_2010 - NYC2$Pop_under_5_CEN_2010 - NYC2$Pop_5_17_CEN_2010;summary(NYC2$over18)
NYC3$over18 <- NYC3$Tot_Population_CEN_2010 - NYC3$Pop_under_5_CEN_2010 - NYC3$Pop_5_17_CEN_2010;summary(NYC3$over18)

#Calls per capita (adults only) - before cutting low population areas. Look at what outliers happen.
NYC1$census_3monthsbefore_pc <- NYC1$census_3monthsbefore/NYC1$over18

NYC1$over18[NYC1$over18==0] <- NA #Make zero population NA
NYC1$census_3monthsbefore_pc <- NYC1$census_3monthsbefore/NYC1$over18

summary(NYC1$census_3monthsbefore_pc) #Max 45
plot(NYC1$census_3monthsbefore_pc)
sd(NYC1$census_3monthsbefore_pc, na.rm=T) #sd 1.29

#Also cut tracts with population below 250 to reduce outliers (as in similar papers).
#Look at the outliers
sum(NYC1$over18<250, na.rm=T)
NYC1[NYC1$over18<250, c("TRACTID","census_3monthsbefore_pc")]
NYC1[NYC1$TRACTID==36061031100,] #Census tract 311 has highest per capita estimate - i I map this in the Census FactFinder, I get a totally normal looking city block, but for some reason in our data the reported population is 2. 
#the 44.5 average is in 36061031100
NYC1$over18[NYC1$over18<250] <- NA
sum(NYC2$over18<250, na.rm=T)
NYC2$over18[NYC2$over18<250] <- NA
sum(NYC3$over18<250, na.rm=T)
NYC3$over18[NYC3$over18<250] <- NA


#Calls per capita (adults only)
NYC1$census_1montharound_pc <- NYC1$census_1montharound/NYC1$over18
#NYC1$census_6monthsbefore_pc <- NYC1$census_6monthsbefore/NYC1$over18 #this won't work as we don't have 6 months out from census (calls start 2010)
NYC1$census_3monthsbefore_pc <- NYC1$census_3monthsbefore/NYC1$over18

NYC2$census_1montharound_pc <- NYC2$census_1montharound/NYC2$over18
#NYC2$census_6monthsbefore_pc <- NYC2$census_6monthsbefore/NYC2$over18
NYC2$census_3monthsbefore_pc <- NYC2$census_3monthsbefore/NYC2$over18

NYC3$census_1montharound_pc <- NYC3$census_1montharound/NYC3$over18
#NYC3$census_6monthsbefore_pc <- NYC3$census_6monthsbefore/NYC3$over18
NYC3$census_3monthsbefore_pc <- NYC3$census_3monthsbefore/NYC3$over18

NYC1$election2010_6monthsbefore_pc <- NYC1$election2010_6monthsbefore/NYC1$over18
NYC2$election2010_6monthsbefore_pc <- NYC2$election2010_6monthsbefore/NYC2$over18
NYC3$election2010_6monthsbefore_pc <- NYC3$election2010_6monthsbefore/NYC3$over18

NYC1$election2010_3monthsbefore_pc <- NYC1$election2010_3monthsbefore/NYC1$over18
NYC2$election2010_3monthsbefore_pc <- NYC2$election2010_3monthsbefore/NYC2$over18
NYC3$election2010_3monthsbefore_pc <- NYC3$election2010_3monthsbefore/NYC3$over18

NYC1$election2010_1montharound_pc <- NYC1$election2010_1montharound/NYC1$over18
NYC2$election2010_1montharound_pc <- NYC2$election2010_1montharound/NYC2$over18
NYC3$election2010_1montharound_pc <- NYC3$election2010_1montharound/NYC3$over18

NYC1$electioncycle_2010_pc <- NYC1$electioncycle_2010/NYC1$over18
NYC2$electioncycle_2010_pc <- NYC2$electioncycle_2010/NYC2$over18
NYC3$electioncycle_2010_pc <- NYC3$electioncycle_2010/NYC3$over18

NYC1$allyear_2010_pc <- NYC1$allyear_2010/NYC1$over18
NYC2$allyear_2010_pc <- NYC2$allyear_2010/NYC2$over18
NYC3$allyear_2010_pc <- NYC3$allyear_2010/NYC3$over18

NYC1$donations2010_count_pc <- NYC1$donations2010_count/NYC1$over18
NYC2$donations2010_count_pc <- NYC2$donations2010_count/NYC2$over18
NYC3$donations2010_count_pc <- NYC3$donations2010_count/NYC3$over18

NYC1$donations2010_dollartotal_pc <- NYC1$donations2010_dollartotal/NYC1$over18
NYC2$donations2010_dollartotal_pc <- NYC2$donations2010_dollartotal/NYC2$over18
NYC3$donations2010_dollartotal_pc <- NYC3$donations2010_dollartotal/NYC3$over18

NYC1$donations2010_3mospreelection_count_pc <- NYC1$donations2010_3mospreelection_count/NYC1$over18
NYC2$donations2010_3mospreelection_count_pc <- NYC2$donations2010_3mospreelection_count/NYC2$over18
NYC3$donations2010_3mospreelection_count_pc <- NYC3$donations2010_3mospreelection_count/NYC3$over18

NYC1$donations2010_6mospreelection_count_pc <- NYC1$donations2010_6mospreelection_count/NYC1$over18
NYC2$donations2010_6mospreelection_count_pc <- NYC2$donations2010_6mospreelection_count/NYC2$over18
NYC3$donations2010_6mospreelection_count_pc <- NYC3$donations2010_6mospreelection_count/NYC3$over18


#Characterize again, now that low population areas out
summary(NYC1$census_3monthsbefore_pc) #Max 0.56
plot(NYC1$census_3monthsbefore_pc)
sd(NYC1$census_3monthsbefore_pc, na.rm=T) #0.04

#How many calls made in 3 month period in places with more than 250 pop?
sum(NYC1$census_3monthsbefore)
sum(NYC1[NYC1$over18>0,"census_3monthsbefore"], na.rm=T)
sum(NYC2$census_3monthsbefore)
sum(NYC2[NYC2$over18>0,"census_3monthsbefore"], na.rm=T)
sum(NYC3$census_3monthsbefore)
sum(NYC3[NYC3$over18>0,"census_3monthsbefore"], na.rm=T)


#Plot calls against other variables

#Census return rates -all
pdf("plot_censusreturn_tracts.pdf")
par(mfrow=c(1,3)) 
plot(NYC1$Mail_Return_Rate_CEN_2010,NYC1$census_1montharound_pc, xlim=c(40,100), xlab="Census Mail Return Rate", ylab="311 calls per capita, 1 month around Census")
plot(NYC1$Mail_Return_Rate_CEN_2010, NYC1$census_3monthsbefore_pc,xlim=c(40,100), xlab="Census Mail Return Rate", ylab="311 calls per capita, 3 months before Census")
dev.off()
#Census return rates - public
pdf("plot_censusreturn_tracts_public.pdf")
par(mfrow=c(1,3)) 
plot(NYC2$Mail_Return_Rate_CEN_2010,NYC2$census_1montharound_pc, xlim=c(40,100), xlab="Census Mail Return Rate", ylab="311 calls per capita, 1 month around Census")
plot(NYC2$Mail_Return_Rate_CEN_2010, NYC2$census_3monthsbefore_pc,xlim=c(40,100), xlab="Census Mail Return Rate", ylab="311 calls per capita, 3 months before Census")
dev.off()
#Census return rates - street
pdf("plot_censusreturn_tracts_street.pdf")
par(mfrow=c(1,3)) 
plot(NYC3$Mail_Return_Rate_CEN_2010,NYC3$census_1montharound_pc, xlim=c(40,100), xlab="Census Mail Return Rate", ylab="311 calls per capita, 1 month around Census")
plot(NYC3$Mail_Return_Rate_CEN_2010, NYC3$census_3monthsbefore_pc,xlim=c(40,100), xlab="Census Mail Return Rate", ylab="311 calls per capita, 3 months before Census")
dev.off()

#Census return rates - all three types, at just 3-month window.
#Checking maximums. They are all in the plot area.
range(NYC1$census_3monthsbefore_pc, na.rm=T) #0.56
range(NYC2$census_3monthsbefore_pc, na.rm=T) #0.47
range(NYC3$census_3monthsbefore_pc, na.rm=T) #0.24

pdf("plot_censusreturn_tracts_allcalltypes.pdf", width=8, height=6)
par(mfrow=c(1,3)) 
plot(NYC1$Mail_Return_Rate_CEN_2010, NYC1$census_3monthsbefore_pc,xlim=c(40,100),ylim=c(0,.6), xlab="Census Mail Return Rate", ylab="All 311 calls per capita, 3 months before Census", main="All 311 calls")
plot(NYC2$Mail_Return_Rate_CEN_2010, NYC2$census_3monthsbefore_pc,xlim=c(40,100),ylim=c(0,.6), xlab="Census Mail Return Rate", ylab="'Public' 311 calls per capita, 3 months before Census", main="'Public' 311 calls")
plot(NYC3$Mail_Return_Rate_CEN_2010, NYC3$census_3monthsbefore_pc,xlim=c(40,100),ylim=c(0,.6), xlab="Census Mail Return Rate", ylab="'Street-level' 311 calls per capita, 3 months before Census", main = "'Street' 311 calls")
dev.off()

pdf("plot_censusreturn_tracts_allcalltypes_addloess.pdf", width=8, height=6)
par(mfrow=c(1,3)) 
library(stats)

scatter.smooth(NYC1$Mail_Return_Rate_CEN_2010, NYC1$census_3monthsbefore_pc,xlim=c(40,100),ylim=c(0,.6), xlab="Census Mail Return Rate", ylab="All 311 calls per capita, 3 months before Census", main="All 311 calls")
scatter.smooth(NYC2$Mail_Return_Rate_CEN_2010, NYC2$census_3monthsbefore_pc,xlim=c(40,100),ylim=c(0,.6), xlab="Census Mail Return Rate", ylab="'Public' 311 calls per capita, 3 months before Census", main="'Public' 311 calls")
scatter.smooth(NYC3$Mail_Return_Rate_CEN_2010, NYC3$census_3monthsbefore_pc,xlim=c(40,100),ylim=c(0,.6), xlab="Census Mail Return Rate", ylab="'Street-level' 311 calls per capita, 3 months before Census", main = "'Street' 311 calls")
dev.off()

plot(NYC1$Mail_Return_Rate_CEN_2010, NYC1$census_3monthsbefore_pc,xlim=c(40,100), xlab="Census Mail Return Rate", ylab="All 311 calls per capita, 3 months before Census", main="All 311 calls")
plot(NYC2$Mail_Return_Rate_CEN_2010, NYC2$census_3monthsbefore_pc,xlim=c(40,100), xlab="Census Mail Return Rate", ylab="'Public' 311 calls per capita, 3 months before Census", main="'Public' 311 calls")
plot(NYC3$Mail_Return_Rate_CEN_2010, NYC3$census_3monthsbefore_pc,xlim=c(40,100), xlab="Census Mail Return Rate", ylab="'Street-level' 311 calls per capita, 3 months before Census", main = "'Street' 311 calls")


#Correlations -all calls
cor(NYC1$census_1montharound_pc, NYC1$Mail_Return_Rate_CEN_2010, use="complete.obs")
census_1month <- lm(NYC1$census_1montharound_pc ~ NYC1$Mail_Return_Rate_CEN_2010)
summary(census_1month)

cor(NYC1$census_3monthsbefore_pc, NYC1$Mail_Return_Rate_CEN_2010, use="complete.obs")
census_3month <- lm(NYC1$census_3monthsbefore_pc ~ NYC1$Mail_Return_Rate_CEN_2010)
summary(census_3month)

#Correlations -public calls
cor(NYC2$census_1montharound_pc, NYC2$Mail_Return_Rate_CEN_2010, use="complete.obs")
census_1month <- lm(NYC2$census_1montharound_pc ~ NYC2$Mail_Return_Rate_CEN_2010)
summary(census_1month)

cor(NYC2$census_3monthsbefore_pc, NYC2$Mail_Return_Rate_CEN_2010, use="complete.obs")
census_3month <- lm(NYC2$census_3monthsbefore_pc ~ NYC2$Mail_Return_Rate_CEN_2010)
summary(census_3month)

#Correlations -street calls
cor(NYC3$census_1montharound_pc, NYC3$Mail_Return_Rate_CEN_2010, use="complete.obs")
census_1month <- lm(NYC3$census_1montharound_pc ~ NYC3$Mail_Return_Rate_CEN_2010)
summary(census_1month)
#instead of small negative correlation, now it's v. v. small negative correlation (basically nothing).
cor(NYC3$census_3monthsbefore_pc, NYC3$Mail_Return_Rate_CEN_2010, use="complete.obs")
census_3month <- lm(NYC3$census_3monthsbefore_pc ~ NYC3$Mail_Return_Rate_CEN_2010)
summary(census_3month)

## look at the same thing but not percapita.
cor(NYC1$census_3monthsbefore, NYC1$Mail_Return_Rate_CEN_2010, use="complete.obs")
cor(NYC2$census_3monthsbefore, NYC2$Mail_Return_Rate_CEN_2010, use="complete.obs")
cor(NYC3$census_3monthsbefore, NYC3$Mail_Return_Rate_CEN_2010, use="complete.obs")

summary(NYC1$College_ACS_08_12)
summary(NYC1$pct_NH_Blk_alone_CEN_2010)
summary(NYC1$Pov_Univ_ACS_08_12)
summary(NYC1$pct_Born_foreign_ACS_08_12)

#output subsets of NYC1,2,3 for plotting
keep <- c("TRACTID", "census_3monthsbefore", "Mail_Return_Rate_CEN_2010", "Tot_Population_CEN_2010", "allyear_2010")

library(foreign)
write.csv(subset(NYC1, select=keep) , file= "tractvolumes_allcalls.csv")

#okay, now rerun these with census vars included.
#first recode some census vars

NYC1$pctpovertyACS <- NYC1$Prs_Blw_Pov_Lev_ACS_08_12/NYC1$Pov_Univ_ACS_08_12
NYC1$pctcollegeACS <- NYC1$College_ACS_08_12/NYC1$Tot_Population_CEN_2010
NYC2$pctpovertyACS <- NYC2$Prs_Blw_Pov_Lev_ACS_08_12/NYC2$Pov_Univ_ACS_08_12
NYC2$pctcollegeACS <- NYC2$College_ACS_08_12/NYC2$Tot_Population_CEN_2010
NYC3$pctpovertyACS <- NYC3$Prs_Blw_Pov_Lev_ACS_08_12/NYC3$Pov_Univ_ACS_08_12
NYC3$pctcollegeACS <- NYC3$College_ACS_08_12/NYC3$Tot_Population_CEN_2010

#first just look at need/calls-- not sure what time window to use.
# start with percent poverty, and look across different timeframes/call type subsets
d1census3mospov <- cor(NYC1$census_3monthsbefore_pc, NYC1$pctpovertyACS, use="complete.obs") 
d1census1mospov <- cor(NYC1$census_1montharound_pc, NYC1$pctpovertyACS, use="complete.obs") 
#d1elec2012pov <- cor(NYC1$electioncycle_2012_pc, NYC1$pctpovertyACS, use="complete.obs") 
d1elect6mospov <- cor(NYC1$election2010_6monthsbefore_pc, NYC1$pctpovertyACS, use="complete.obs") 
d1elect1mospov <- cor(NYC1$election2010_1montharound_pc, NYC1$pctpovertyACS, use="complete.obs") 
d1elect3mospov <- cor(NYC1$election2010_3monthsbefore_pc, NYC1$pctpovertyACS, use="complete.obs") 

d2census3mospov <-cor(NYC2$census_3monthsbefore_pc, NYC2$pctpovertyACS, use="complete.obs") 
d2census1mospov <-cor(NYC2$census_1montharound_pc, NYC2$pctpovertyACS, use="complete.obs") 
#d2elec2012pov <- cor(NYC2$electioncycle_2012_pc, NYC2$pctpovertyACS, use="complete.obs") 
d2elect6mospov <-cor(NYC2$election2010_6monthsbefore_pc, NYC2$pctpovertyACS, use="complete.obs") 
d2elect1mospov <- cor(NYC2$election2010_1montharound_pc, NYC2$pctpovertyACS, use="complete.obs") 
d2elect3mospov <- cor(NYC2$election2010_3monthsbefore_pc, NYC2$pctpovertyACS, use="complete.obs") 

d3census3mospov <-cor(NYC3$census_3monthsbefore_pc, NYC3$pctpovertyACS, use="complete.obs") 
d3census1mospov <-cor(NYC3$census_1montharound_pc, NYC3$pctpovertyACS, use="complete.obs") 
#d3elec2012pov <- cor(NYC3$electioncycle_2012_pc, NYC3$pctpovertyACS, use="complete.obs") 
d3elect6mospov <-cor(NYC3$election2010_6monthsbefore_pc, NYC3$pctpovertyACS, use="complete.obs") 
d3elect1mospov <- cor(NYC3$election2010_1montharound_pc, NYC3$pctpovertyACS, use="complete.obs") 
d3elect3mospov <- cor(NYC3$election2010_3monthsbefore_pc, NYC3$pctpovertyACS, use="complete.obs") 

#now switch to median HH income, again looking at different times/calls.
NYC1$medhhinc <- as.numeric(gsub("[,]", "", gsub("[$]", "", as.character(NYC1$Med_HHD_Inc_ACS_08_12)))) 
NYC2$medhhinc <- as.numeric(gsub("[,]", "", gsub("[$]", "", as.character(NYC2$Med_HHD_Inc_ACS_08_12))))
NYC3$medhhinc <- as.numeric(gsub("[,]", "", gsub("[$]", "", as.character(NYC3$Med_HHD_Inc_ACS_08_12))))

d1census3mosinc <- cor(NYC1$census_3monthsbefore_pc, NYC1$medhhinc, use="complete.obs") 
d1census1mosinc <- cor(NYC1$census_1montharound_pc, NYC1$medhhinc, use="complete.obs")
d1elect6mosinc <- cor(NYC1$election2010_6monthsbefore_pc, NYC1$medhhinc, use="complete.obs") 
d1elect1mosinc <- cor(NYC1$election2010_1montharound_pc, NYC1$medhhinc, use="complete.obs") 
d1elect3mosinc <- cor(NYC1$election2010_3monthsbefore_pc, NYC1$medhhinc, use="complete.obs") 

d2census3mosinc <-cor(NYC2$census_3monthsbefore_pc, NYC2$medhhinc, use="complete.obs") 
d2census1mosinc  <-cor(NYC2$census_1montharound_pc, NYC2$medhhinc, use="complete.obs") 
d2elect6mosinc <-cor(NYC2$election2010_6monthsbefore_pc, NYC2$medhhinc, use="complete.obs") 
d2elect1mosinc <- cor(NYC2$election2010_1montharound_pc, NYC2$medhhinc, use="complete.obs") 
d2elect3mosinc <- cor(NYC2$election2010_3monthsbefore_pc, NYC2$medhhinc, use="complete.obs") 

d3census3mosinc  <-cor(NYC3$census_3monthsbefore_pc, NYC3$medhhinc, use="complete.obs") 
d3census1mosinc <-cor(NYC3$census_1montharound_pc, NYC3$medhhinc, use="complete.obs") 
d3elect6mosinc  <-cor(NYC3$election2010_6monthsbefore_pc, NYC3$medhhinc, use="complete.obs") 
d3elect1mosinc  <- cor(NYC3$election2010_1montharound_pc, NYC3$medhhinc, use="complete.obs") 
d3elect3mosinc  <- cor(NYC3$election2010_3monthsbefore_pc, NYC3$medhhinc, use="complete.obs") 

### so: not too reliable.  

## try plotting this. 
pdf(file="needand311_poverty_alltimesallcalls.pdf")
plot(NULL, NULL, ylim=c(-.2, .3), xlim=c(.5,3.5), main="Tract-level Poverty and 311 Calls Per Capita", ylab= "Correlation", xlab="Timeframe", xaxt="n", cex.axis=1.2, cex.lab=1.2, cex.main=1.5)
axis(side = 1, at = c(1,2,3), labels = c("1 Month", "3 Months", "6 Months"), cex.axis=1.2)
abline(h=0, lty=2, col=gray(.5))
points(1, d1elect1mospov, pch=1, cex=2, lwd=2)
points(1, d1elect3mospov, pch=2, cex=2, lwd=2)
points(1, d1elect6mospov, pch=0, cex=2, lwd=2)
points(2, d2elect1mospov, pch=1, cex=2, lwd=2)
points(2, d2elect3mospov, pch=2, cex=2, lwd=2)
points(2, d2elect6mospov, pch=0, cex=2, lwd=2)
points(3, d3elect1mospov, pch=1, cex=2, lwd=2)
points(3, d3elect3mospov, pch=2, cex=2, lwd=2)
points(3, d3elect6mospov, pch=0, cex=2, lwd=2)
legend("topright", legend = c("All Calls", "`Public'", "`Street-level'"), pch=c(1,2,0), lwd=2, lty=NA, cex=1.5)
dev.off()

pdf(file="needand311_medHHincome_alltimesallcalls.pdf")
plot(NULL, NULL, ylim=c(-.2, .3), xlim=c(.5,3.5), main="Tract-level Median Income and 311 Calls Per Capita", ylab= "Correlation", xlab="Timeframe", xaxt="n", cex.axis=1.2, cex.lab=1.2, cex.main=1.5)
axis(side = 1, at = c(1,2,3), labels = c("1 Month", "3 Months", "6 Months"), cex.axis=1.2)
abline(h=0, lty=2, col=gray(.5))
points(1, d1elect1mosinc, pch=1, cex=2, lwd=2)
points(1, d1elect3mosinc, pch=2, cex=2, lwd=2)
points(1, d1elect6mosinc, pch=0, cex=2, lwd=2)
points(2, d2elect1mosinc, pch=1, cex=2, lwd=2)
points(2, d2elect3mosinc, pch=2, cex=2, lwd=2)
points(2, d2elect6mosinc, pch=0, cex=2, lwd=2)
points(3, d3elect1mosinc, pch=1, cex=2, lwd=2)
points(3, d3elect3mosinc, pch=2, cex=2, lwd=2)
points(3, d3elect6mosinc, pch=0, cex=2, lwd=2)
legend("bottomright", legend = c("All Calls", "`Public'", "`Street-level'"), pch=c(1,2,0), lwd=2, lty=NA, cex=1.5)
dev.off()


census_1month <- lm(NYC1$census_1montharound_pc ~ NYC1$Mail_Return_Rate_CEN_2010)
summary(census_1month)
census_1monthcov <- lm(census_1montharound_pc ~ Mail_Return_Rate_CEN_2010+ pctpovertyACS + pctcollegeACS + pct_NH_Blk_alone_CEN_2010 + pct_Born_foreign_ACS_08_12, data=NYC1)
summary(census_1monthcov) #doesn't change much at all.

census_3month_1 <- lm(census_3monthsbefore_pc ~ Mail_Return_Rate_CEN_2010, data=NYC1)
summary(census_3month_1)
census_3monthneed_1 <- lm(census_3monthsbefore_pc ~ Mail_Return_Rate_CEN_2010+pctpovertyACS +medhhinc, data=NYC1)
summary(census_3monthneed_1) 
census_3monthcov_1 <- lm(census_3monthsbefore_pc ~ Mail_Return_Rate_CEN_2010+ pctpovertyACS +medhhinc+ pctcollegeACS + pct_NH_Blk_alone_CEN_2010 + pct_Born_foreign_ACS_08_12, data=NYC1)
summary(census_3monthcov_1)

#other subsets
census_3month_2 <- lm(census_3monthsbefore_pc ~ Mail_Return_Rate_CEN_2010, data=NYC2)
summary(census_3month_2)
census_3monthneed_2 <- lm(census_3monthsbefore_pc ~ Mail_Return_Rate_CEN_2010+ pctpovertyACS +medhhinc, data=NYC2)
summary(census_3monthneed_2) 
census_3monthcov_2 <- lm(census_3monthsbefore_pc ~ Mail_Return_Rate_CEN_2010+ pctpovertyACS +medhhinc+ pctcollegeACS + pct_NH_Blk_alone_CEN_2010 + pct_Born_foreign_ACS_08_12, data=NYC2)
summary(census_3monthcov_2) #

census_3month_3 <- lm(census_3monthsbefore_pc ~ Mail_Return_Rate_CEN_2010, data=NYC3)
summary(census_3month_3)
census_3monthneed_3 <- lm(census_3monthsbefore_pc ~ Mail_Return_Rate_CEN_2010+ pctpovertyACS +medhhinc, data=NYC3)
summary(census_3monthneed_3) 
census_3monthcov_3 <- lm(census_3monthsbefore_pc ~ Mail_Return_Rate_CEN_2010+ pctpovertyACS +medhhinc+ pctcollegeACS + pct_NH_Blk_alone_CEN_2010 + pct_Born_foreign_ACS_08_12, data=NYC3)
summary(census_3monthcov_3) 

#adding in pct black (census), pct in poverty (ACS), pct foreign born (ACS), pct college degree (ACS) doesn't change the fundamental (non) relationship between census returns and 311 

## now, make one panel of what will ultimately be a three-panel plot, demonstrating how adding in covars doesn't really matter.
allcallsest <- c(coef(census_3month_1)[2],coef(census_3monthneed_1)[2],coef(census_3monthcov_1)[2])
censuspanel_allcalls <- as.data.frame(cbind(allcallsest, rbind(confint(census_3month_1, 'Mail_Return_Rate_CEN_2010', level=0.95), 
confint(census_3monthneed_1, 'Mail_Return_Rate_CEN_2010', level=0.95),
confint(census_3monthcov_1, 'Mail_Return_Rate_CEN_2010', level=0.95))))
colnames(censuspanel_allcalls) <- c("est", "CIlow", "CIhigh")

pdf("addingcovars_censusresponse_allcalls_3mos.pdf")
plot(c(1:3), censuspanel_allcalls$est, xlim=c(.5, 3.5), ylim=c(-.01, .01), ylab="Coefficient on Census Response Rate", xlab="", xaxt='n', main="Census Response and 311 Calls", pch=16, cex=1.1)
abline(h=0, lty=2, col=gray(.5))
axis(side = 1, at = c(1,2,3), labels = c("No Covars", "Need (Poverty, HH Income)", "Add'l Covariates"))
segments(c(1:3), censuspanel_allcalls$CIlow, c(1:3),censuspanel_allcalls$CIhigh, lwd=1.5)
dev.off()

#output tables for these models as well.
library(stargazer)
stargazer(census_3month_1,census_3monthneed_1,census_3monthcov_1, out="addingcovars_censusresponse_allcalls_3mos.tex", 
	title = "Regressing Tract-level 311 Call Volumes onto Census Response Rates and Covariates",
          covariate.labels = c("Census Return Rate", "Percent Poverty (ACS)", "Median HH Income", "Percent College (ACS)", 
"Percent Black", "Percent Foreign Born (ACS)"),
          dep.var.labels   = "311 Calls Per Capita", label="censuscovars", omit.stat = "f") 


#replicate this using raw #s and a pop. control instead of per capita (per reviewer's question)

census_3month_1 <- lm(census_3monthsbefore ~ Mail_Return_Rate_CEN_2010 + over18, data=NYC1)
summary(census_3month_1)
census_3monthneed_1 <- lm(census_3monthsbefore ~ Mail_Return_Rate_CEN_2010 + over18 +pctpovertyACS +medhhinc, data=NYC1)
summary(census_3monthneed_1) 
census_3monthcov_1 <- lm(census_3monthsbefore ~ Mail_Return_Rate_CEN_2010+ over18 + pctpovertyACS +medhhinc+ pctcollegeACS + pct_NH_Blk_alone_CEN_2010 + pct_Born_foreign_ACS_08_12, data=NYC1)
summary(census_3monthcov_1)

library(stargazer)
stargazer(census_3month_1,census_3monthneed_1,census_3monthcov_1, out="REPLICATION_addingcovars_censusresponse_allcalls_3mos.tex", 
	title = "Regressing Tract-level 311 Call Volumes onto Census Response Rates and Covariates (Alternative Population Approach)",
          covariate.labels = c("Census Return Rate", "Over-18 Population", "Percent Poverty (ACS)", "Median HH Income", "Percent College (ACS)", 
"Percent Black", "Percent Foreign Born (ACS)"),
          dep.var.labels   = "311 Calls Per Capita", label="rep_censuscovars", omit.stat = "f") 


############################
#okay, now switching over to look at political donations, in about the same ways

## make a combined plot: looking at donation counts for some time frame, per capita, across all three types of calls in a panel (like the other plots in the paper)

pdf("plot_donationcounts_tracts_allcalltypes.pdf", width=8, height=6)
par(mfrow=c(1,3)) 

plot(NYC1$donations2010_count_pc, NYC1$electioncycle_2010_pc, ylim=c(0,1.5), xlab="Individual Political Donations Reported per capita", ylab="311 calls per capita, 2010 Election Cycle (2009-2010)", main="All 311 calls")
plot(NYC2$donations2010_count_pc, NYC2$electioncycle_2010_pc, ylim=c(0,1.5), xlab="Individual Political Donations Reported per capita", ylab="`Public` 311 calls per capita, 2010 Election Cycle (2009-2010)", main="'Public' 311 calls")
plot(NYC3$donations2010_count_pc, NYC3$electioncycle_2010_pc, ylim=c(0,1.5), xlab="Individual Political Donations Reported per capita", ylab="`Street-level` 311 calls per capita, 2010 Election Cycle (2009-2010)", main = "'Street' 311 calls")

dev.off()

pdf("plot_donationcounts_tracts_allcalltypes_3months.pdf", width=8, height=6)
par(mfrow=c(1,3)) 

plot(NYC1$donations2010_3mospreelection_count_pc, NYC1$electioncycle_2010_pc, ylim=c(0,1.5), xlab="Individual Political Donations Reported per capita", ylab="311 calls per capita, 3 months preceding 2010 election", main="All 311 calls")
plot(NYC2$donations2010_3mospreelection_count_pc, NYC2$electioncycle_2010_pc, ylim=c(0,1.5), xlab="Individual Political Donations Reported per capita", ylab="`Public` 311 calls per capita, 3 months preceding 2010 election", main="'Public' 311 calls")
plot(NYC3$donations2010_3mospreelection_count_pc, NYC3$electioncycle_2010_pc, ylim=c(0,1.5), xlab="Individual Political Donations Reported per capita", ylab="`Street-level` 311 calls per capita, 3 months preceding 2010 election", main = "'Street' 311 calls")
dev.off()


pdf("plot_donationcounts_tracts_allcalltypes_6months.pdf", width=8, height=6)
par(mfrow=c(1,3)) 

plot(NYC1$donations2010_6mospreelection_count_pc, NYC1$election2010_6monthsbefore_pc, ylim=c(0,1.5), xlab="Individual Political Donations Reported per capita", ylab="311 calls per capita, 6 months preceding 2010 election", main="All 311 calls")
plot(NYC2$donations2010_6mospreelection_count_pc, NYC2$election2010_6monthsbefore_pc, ylim=c(0,1.5), xlab="Individual Political Donations Reported per capita", ylab="`Public` 311 calls per capita, 6 months preceding 2010 election", main="'Public' 311 calls")
plot(NYC3$donations2010_6mospreelection_count_pc, NYC3$election2010_6monthsbefore_pc, ylim=c(0,1.5), xlab="Individual Political Donations Reported per capita", ylab="`Street-level` 311 calls per capita, 6 months preceding 2010 election", main = "'Street' 311 calls")
dev.off()


#Correlations -all calls
cor(NYC1$donations2010_3mospreelection_count, NYC1$electioncycle_2010_pc, use="complete.obs")
donations_fullcycle <- lm(NYC1$electioncycle_2010_pc ~ NYC1$donations2010_3mospreelection_count)
summary(donations_fullcycle)

	#per capita
cor(NYC1$donations2010_3mospreelection_count_pc, NYC1$electioncycle_2010_pc, use="complete.obs")
donations_fullcycle <- lm(NYC1$electioncycle_2010_pc ~ NYC1$donations2010_3mospreelection_count_pc)
summary(donations_fullcycle)

donationspc_dropoutliers <- NYC1[NYC1$donations2010_3mospreelection_count_pc < .065,] #pretty much eyeballing this from the plot
dim(NYC1); dim(donationspc_dropoutliers)
donations_fullcycle_dropoutliers <- lm(electioncycle_2010_pc ~ donations2010_3mospreelection_count_pc, data= donationspc_dropoutliers )
summary(donations_fullcycle_dropoutliers)
plot(donationspc_dropoutliers$donations2010_3mospreelection_count_pc, donationspc_dropoutliers$electioncycle_2010_pc, ylim=c(0,3), xlab="Individual Political Donations Reported per capita", ylab="311 calls per capita, 2010 Election Cycle (2011-2010)")
cor(donationspc_dropoutliers$donations2010_3mospreelection_count_pc, donationspc_dropoutliers$electioncycle_2010_pc, use="complete.obs")

	#per capita including 0
NYC1$donations2010_3mospreelection_count_pc0 <- NYC1$donations2010_3mospreelection_count_pc
NYC1[is.na(NYC1$donations2010_3mospreelection_count_pc), "donations2010_3mospreelection_count_pc0"] <- 0
summary(NYC1$donations2010_3mospreelection_count_pc); summary(NYC1$donations2010_3mospreelection_count_pc0)
cor(NYC1$donations2010_3mospreelection_count_pc, NYC1$electioncycle_2010_pc, use="complete.obs")
cor(NYC1$donations2010_3mospreelection_count_pc0, NYC1$electioncycle_2010_pc, use="complete.obs")
donations_fullcycle <- lm(NYC1$electioncycle_2010_pc ~ NYC1$donations2010_3mospreelection_count_pc)
summary(donations_fullcycle)


	#dollar value per capita
cor(NYC1$donations2010_dollartotal_pc, NYC1$electioncycle_2010_pc, use="complete.obs")
donationdollars_fullcycle <- lm(NYC1$electioncycle_2010_pc ~ NYC1$donations2010_dollartotal_pc)
summary(donationdollars_fullcycle)

donationdollarspc_dropoutliers <- NYC1[NYC1$donations2010_dollartotal_pc < 500,] #pretty much eyeballing this from the plot
dim(NYC1); dim(donationdollarspc_dropoutliers)
donationdollars_fullcycle_dropoutliers <- lm(electioncycle_2010_pc ~ donations2010_dollartotal_pc, data= donationdollarspc_dropoutliers )
summary(donationdollars_fullcycle_dropoutliers) #less exciting, but still there
plot(donationdollarspc_dropoutliers$donations2010_dollartotal_pc, donationdollarspc_dropoutliers$electioncycle_2010_pc, ylim=c(0,3), xlab="Value of Individual Political Donations Reported per capita", ylab="311 calls per capita, 2010 Election Cycle (2011-2010)")
cor(donationdollarspc_dropoutliers$donations2010_3mospreelection_count_pc, donationdollarspc_dropoutliers$electioncycle_2010_pc, use="complete.obs")


#Correlations -public calls
	#donations per capita
cor(NYC2$donations2010_3mospreelection_count_pc, NYC2$electioncycle_2010_pc, use="complete.obs")
donations_fullcycle_public <- lm(NYC2$electioncycle_2010_pc ~ NYC2$donations2010_3mospreelection_count_pc)
summary(donations_fullcycle_public)
	#dollar value per capita
cor(NYC2$donations2010_dollartotal_pc, NYC2$electioncycle_2010_pc, use="complete.obs")
donationdollars_fullcycle_public <- lm(NYC2$electioncycle_2010_pc ~ NYC2$donations2010_dollartotal_pc)
summary(donationdollars_fullcycle_public)


#Correlations -street calls
	#donations per capita
cor(NYC3$donations2010_3mospreelection_count_pc, NYC3$electioncycle_2010_pc, use="complete.obs")
donations_fullcycle_street <- lm(NYC3$electioncycle_2010_pc ~ NYC3$donations2010_3mospreelection_count_pc)
summary(donations_fullcycle_street)

	#dollar value per capita
cor(NYC3$donations2010_dollartotal_pc, NYC3$electioncycle_2010_pc, use="complete.obs")
donationdollars_fullcycle_street <- lm(NYC3$electioncycle_2010_pc ~ NYC3$donations2010_dollartotal_pc)
summary(donationdollars_fullcycle_street)


###################################################
# down here, make a combined table: all three outcomes with and w/o covars.
#actually, make three of these (for the three subsets of the data), but start with street-level calls.

NYC1$medhhinc <- as.numeric(gsub("[,]", "", gsub("[$]", "", as.character(NYC1$Med_HHD_Inc_ACS_08_12)))) #also want to include this
NYC2$medhhinc <- as.numeric(gsub("[,]", "", gsub("[$]", "", as.character(NYC2$Med_HHD_Inc_ACS_08_12))))
NYC3$medhhinc <- as.numeric(gsub("[,]", "", gsub("[$]", "", as.character(NYC3$Med_HHD_Inc_ACS_08_12))))

donations_3month_1 <- lm(election2010_3monthsbefore_pc ~ donations2010_3mospreelection_count_pc, data=NYC1)
summary(donations_3month_1)
donations_3month_need_1 <- lm(election2010_3monthsbefore_pc ~ donations2010_3mospreelection_count_pc+ pctpovertyACS + medhhinc , data=NYC1)
summary(donations_3month_need_1)
donations_3monthcov_1 <- lm(election2010_3monthsbefore_pc ~ donations2010_3mospreelection_count_pc+ pctpovertyACS + medhhinc + pctcollegeACS + pct_NH_Blk_alone_CEN_2010 + pct_Born_foreign_ACS_08_12, data=NYC1)
summary(donations_3monthcov_1) 

donations_3month_2 <- lm(election2010_3monthsbefore_pc ~ donations2010_3mospreelection_count_pc, data=NYC2)
summary(donations_3month_2)
donations_3month_need_2 <- lm(election2010_3monthsbefore_pc ~ donations2010_3mospreelection_count_pc+ pctpovertyACS + medhhinc , data=NYC2)
summary(donations_3month_need_2)
donations_3monthcov_2 <- lm(election2010_3monthsbefore_pc ~ donations2010_3mospreelection_count_pc+ pctpovertyACS + medhhinc + pctcollegeACS + pct_NH_Blk_alone_CEN_2010 + pct_Born_foreign_ACS_08_12, data=NYC2)
summary(donations_3monthcov_2) 

donations_3month_3 <- lm(election2010_3monthsbefore_pc ~ donations2010_3mospreelection_count_pc, data=NYC3)
summary(donations_3month_3)
donations_3month_need_3 <- lm(election2010_3monthsbefore_pc ~ donations2010_3mospreelection_count_pc+ pctpovertyACS + medhhinc , data=NYC3)
summary(donations_3month_need_3)
donations_3monthcov_3 <- lm(election2010_3monthsbefore_pc ~ donations2010_3mospreelection_count_pc+ pctpovertyACS + medhhinc + pctcollegeACS + pct_NH_Blk_alone_CEN_2010 + pct_Born_foreign_ACS_08_12, data=NYC3)
summary(donations_3monthcov_3) 

#and do this first one with 0s included as well
donations_3month_1_0 <- lm(election2010_3monthsbefore_pc ~ donations2010_3mospreelection_count_pc0, data=NYC1)
summary(donations_3month_1_0)
donations_3month_need_1_0 <- lm(election2010_3monthsbefore_pc ~ donations2010_3mospreelection_count_pc0+ pctpovertyACS + medhhinc , data=NYC1)
summary(donations_3month_need_1_0)
donations_3monthcov_1_0 <- lm(election2010_3monthsbefore_pc ~ donations2010_3mospreelection_count_pc0+ pctpovertyACS + medhhinc + pctcollegeACS + pct_NH_Blk_alone_CEN_2010 + pct_Born_foreign_ACS_08_12, data=NYC1)
summary(donations_3monthcov_1_0) 

## now, make one panel of what will ultimately be a three-panel plot, demonstrating how adding in covars doesn't really matter.
donationsallcallsest <- c(coef(donations_3month_1)[2],coef(donations_3month_need_1)[2],coef(donations_3monthcov_1)[2])
donations_allcalls <- as.data.frame(cbind(donationsallcallsest, rbind(confint(donations_3month_1, 'donations2010_3mospreelection_count_pc', level=0.95), 
confint(donations_3month_need_1, 'donations2010_3mospreelection_count_pc', level=0.95),
confint(donations_3monthcov_1, 'donations2010_3mospreelection_count_pc', level=0.95))))
colnames(donations_allcalls) <- c("est", "CIlow", "CIhigh"); rownames(donations_allcalls)<- c("raw", "need", "covar")

pdf("addingcovars_poldonations_allcalls_3mos.pdf")
plot(c(1:3), donations_allcalls$est, xlim=c(.5, 3.5), ylim=c(-.1, 2), ylab="Coefficient on Total Donations Per Capita (3 months preceding 2010 election)", xlab="", xaxt='n', main="Political Donations and 311 Calls", pch=16, cex=1.1)
abline(h=0, lty=2, col=gray(.5))
axis(side = 1, at = c(1,2,3), labels = c("No Covars", "Need (Poverty, HH Income)", "Add'l Covariates"))
segments(c(1:3), donations_allcalls$CIlow, c(1:3),donations_allcalls$CIhigh, lwd=1.5)
dev.off()


## okay, and make these tables as well.
library(stargazer)
stargazer(donations_3month_1,donations_3month_need_1,donations_3monthcov_1, out="addingcovars_donations_allcalls_3mos.tex", 
	title = "Regressing Tract-level 311 Call Volumes onto Political Donation Rates and Covariates" ,
          covariate.labels = c("Per Capita Donations (Count)", "Percent Poverty (ACS)", "Median HH Income", "Percent College (ACS)", 
"Percent Black", "Percent Foreign Born (ACS)"),
          dep.var.labels   = "311 Calls Per Capita", label="donationcovars", omit.stat = "f") 

######### and also do all this including 0s.
## now, make one panel of what will ultimately be a three-panel plot, demonstrating how adding in covars doesn't really matter.
donationsallcallsest0 <- c(coef(donations_3month_1_0)[2],coef(donations_3month_need_1_0)[2],coef(donations_3monthcov_1_0)[2])
donations_allcalls0 <- as.data.frame(cbind(donationsallcallsest0, rbind(confint(donations_3month_1_0, 'donations2010_3mospreelection_count_pc0', level=0.95), 
confint(donations_3month_need_1_0, 'donations2010_3mospreelection_count_pc0', level=0.95),
confint(donations_3monthcov_1_0, 'donations2010_3mospreelection_count_pc0', level=0.95))))
colnames(donations_allcalls0) <- c("est", "CIlow", "CIhigh"); rownames(donations_allcalls0)<- c("raw", "need", "covar")

pdf("addingcovars_poldonations_allcalls_3mos_with0.pdf")
plot(c(1:3), donations_allcalls0$est, xlim=c(.5, 3.5), ylim=c(-.1, 2), ylab="Coefficient on Total Donations Per Capita (3 months preceding 2010 election)", xlab="", xaxt='n', main="Political Donations and 311 Calls", pch=16, cex=1.1)
abline(h=0, lty=2, col=gray(.5))
axis(side = 1, at = c(1,2,3), labels = c("No Covars", "Need (Poverty, HH Income)", "Add'l Covariates"))
segments(c(1:3), donations_allcalls0$CIlow, c(1:3),donations_allcalls0$CIhigh, lwd=1.5)
dev.off()


## okay, and make these tables as well.
library(stargazer)
stargazer(donations_3month_1_0,donations_3month_need_1_0,donations_3monthcov_1_0, out="addingcovars_donations_allcalls_3mos_with0.tex", 
	title = "Regressing Tract-level 311 Call Volumes onto Political Donation Rates and Covariates" ,
          covariate.labels = c("Per Capita Donations (Count)", "Percent Poverty (ACS)", "Median HH Income", "Percent College (ACS)", 
"Percent Black", "Percent Foreign Born (ACS)"),
          dep.var.labels   = "311 Calls Per Capita", label="donationcovars", omit.stat = "f") 





##now make the replication version of this table for reviewer 3 (with the pop control instead of per capita)
#what's the range of population?
summary(NYC1$over18)

donations_3month_1 <- lm(election2010_3monthsbefore ~ donations2010_3mospreelection_count_pc + over18, data=NYC1)
summary(donations_3month_1)
donations_3month_need_1 <- lm(election2010_3monthsbefore ~ donations2010_3mospreelection_count_pc + over18 + pctpovertyACS + medhhinc, data=NYC1)
summary(donations_3month_need_1)
donations_3monthcov_1 <- lm(election2010_3monthsbefore ~ donations2010_3mospreelection_count_pc + over18 + pctpovertyACS + medhhinc + pctcollegeACS + pct_NH_Blk_alone_CEN_2010 + pct_Born_foreign_ACS_08_12, data=NYC1)
summary(donations_3monthcov_1) 

library(stargazer)
stargazer(donations_3month_1,donations_3month_need_1,donations_3monthcov_1, out="REPLICATION_addingcovars_donations_allcalls_3mos.tex", 
	title = "Regressing Tract-level 311 Call Volumes onto Political Donation Rates and Covariates (Alternative Population Approach)" ,
          covariate.labels = c("Per Capita Donations (Count)", "Over-18 Population", "Percent Poverty (ACS)", "Median HH Income", "Percent College (ACS)", 
"Percent Black", "Percent Foreign Born (ACS)"),
          dep.var.labels   = "311 Calls Per Capita", label="rep_donationcovars", omit.stat = "f") 


#and do it dropping down to raw donation counts as well
donations_3month_1 <- lm(election2010_3monthsbefore ~ donations2010_3mospreelection_count + over18, data=NYC1)
summary(donations_3month_1)
donations_3month_need_1 <- lm(election2010_3monthsbefore ~ donations2010_3mospreelection_count + over18 + Pov_Univ_ACS_08_12 , data=NYC1)
summary(donations_3month_need_1)
donations_3monthcov_1 <- lm(election2010_3monthsbefore ~ donations2010_3mospreelection_count + over18 + Pov_Univ_ACS_08_12  + College_ACS_08_12 + NH_Blk_alone_CEN_2010 + Born_foreign_ACS_08_12, data=NYC1)
summary(donations_3monthcov_1) 

library(stargazer)
stargazer(donations_3month_1,donations_3month_need_1,donations_3monthcov_1, out="REPLICATION1_addingcovars_donations_allcalls_3mos.tex", 
	title = "Regressing Tract-level 311 Call Volumes onto Political Donation Rates and Covariates (Alternative Population Approach)" ,
          covariate.labels = c("Donations (Count)", "Over-18 Population", "Number in Poverty (ACS)", "Number College (ACS)", 
"Number Black", "Number Foreign Born (ACS)"),
          dep.var.labels   = "311 Calls", label="rep1_donationcovars", omit.stat = "f") 

#pull in 0 counts.
NYC1$donations2010_3mospreelection_count0 <- NYC1$donations2010_3mospreelection_count
NYC1[is.na(NYC1$donations2010_3mospreelection_count), "donations2010_3mospreelection_count0"] <- 0
summary(NYC1$donations2010_3mospreelection_count0)
donations_3month_1 <- lm(election2010_3monthsbefore ~ donations2010_3mospreelection_count0 + over18, data=NYC1)
summary(donations_3month_1)
donations_3month_need_1 <- lm(election2010_3monthsbefore ~ donations2010_3mospreelection_count0 + over18 + Pov_Univ_ACS_08_12 , data=NYC1)
summary(donations_3month_need_1)
donations_3monthcov_1 <- lm(election2010_3monthsbefore ~ donations2010_3mospreelection_count0 + over18 + Pov_Univ_ACS_08_12  + College_ACS_08_12 + NH_Blk_alone_CEN_2010 + Born_foreign_ACS_08_12, data=NYC1)
summary(donations_3monthcov_1) 

library(stargazer)
stargazer(donations_3month_1,donations_3month_need_1,donations_3monthcov_1, out="REPLICATION2_addingcovars_donations_allcalls_3mos.tex", 
	title = "Regressing Tract-level 311 Call Volumes onto Political Donation Rates and Covariates (Alternative Population Approach)" ,
          covariate.labels = c("Donations (Count)", "Over-18 Population", "Number in Poverty (ACS)", "Number College (ACS)", 
"Number Black", "Number Foreign Born (ACS)"),
          dep.var.labels   = "311 Calls", label="rep2_donationcovars", omit.stat = "f") 



